A rigorous theoretical and numerical analysis of a nonlinear reaction-diffusion epidemic model pertaining dynamics of COVID-19

The spatial movement of the human population from one region to another and the existence of super-spreaders are the main factors that enhanced the disease incidence. Super-spreaders refer to the individuals having transmitting ability to multiple pathogens. In this article, an epidemic model with spatial and temporal effects is formulated to analyze the impact of some preventing measures of COVID-19. The model is developed using six nonlinear partial differential equations. The infectious individuals are sub-divided into symptomatic, asymptomatic and super-spreader classes. In this study, we focused on the rigorous qualitative analysis of the reaction-diffusion model. The fundamental mathematical properties of the proposed COVID-19 epidemic model such as boundedness, positivity, and invariant region of the problem solution are derived, which ensure the validity of the proposed model. The model equilibria and its stability analysis for both local and global cases have been presented. The normalized sensitivity analysis of the model is carried out in order to observe the crucial factors in the transmission of infection. Furthermore, an efficient numerical scheme is applied to solve the proposed model and detailed simulation are performed. Based on the graphical observation, diffusion in the context of confined public gatherings is observed to significantly inhibit the spread of infection when compared to the absence of diffusion. This is especially important in scenarios where super-spreaders may play a major role in transmission. The impact of some non-pharmaceutical interventions are illustrated graphically with and without diffusion. We believe that the present investigation will be beneficial in understanding the complex dynamics and control of COVID-19 under various non-pharmaceutical interventions.

It has been declared a pandemic by World Health Organization, as it spread very rapidly among the human population through different sources and reported million of confirmed cases accompanied by millions of deaths throughout the world 1 .The transmission of this infectious disease is difficult to control due to the uncertain nature of the virus.Many countries implement social distancing policies and avoid public gatherings, isolating the infected individuals to control the disease incidence.It is still a major threat to public health, although several vaccines are available now.The major factors that rapid the infection transmission are the super-spreaders and spatial movement of populations, since the disease may be transmitted faster in one place than another because of social contacts 2 .In earlier studies of disease epidemiology, it was assumed that susceptible hosts within a population had equal chances to become infected 3 .Further studies uncover the fact that heterogeneities in pathogen • Each newborn can get the infection.The susceptible class increased with newborns and reduced through infection and natural death.• The susceptible get an infection after interacting with infected individuals and are moved to the exposed class which enhances it while reducing natural death and population completing their latency/incubation period.• The fraction of the exposed population that has clinical symptoms is moved to the symptomatically infected class.This class is reduced with natural death, death due to infection, and recovery from infection.• The fraction of the exposed population that has the ability to transfer multiple pathogens are super-spreaders and added to class super-spreaders.This class is also reduced with natural death, death due to infection, and recovery from infection.• The fraction of the exposed population that has no clinical symptoms is moved to the asymptomatically infected class, which reduces with natural death and recovery from infection.• The recovered class varies by moving individuals recovered from infection in any of the respective compart- ments and natural death.
Considering the above listed assumptions the spatio-temporal compartmental model describing the dynamics of COVID-19 is described as follows: where ( t, ỹ) ∈ [0, T max ] × [a, b] ; T max > 0 and denote the force of infection, which represents the transmission potential when susceptible individuals interact with infectious individuals I 1 , I 2 , and I 3 .The detailed description of the parameters in model (1) is tabulated in Table 2, the coefficients of diffusivity are denoted by D i for i = 1, 2, . . ., 6 .Moreover, the following no-flux boundary conditions are considered for problem 1: The symbol ξ denotes the state variables of the model (1), where in (2), ξ ỹ represent partial derivatives of each state variable of model (1) with respect to the spatial variable ỹ.

Initial conditions
To simulate model (1), the initial conditions (ICs) given by ( 3) and ( 4) are used.The ICs (4) are chosen based on 27 : (1) ξ ỹ ( t, −2) = 0, ξ ỹ ( t, 2) = 0. where S 0 = 34, 806, 871, E 0 = 20, 000, I 10 = 3.0,I 20 = 110.0,I 30 = 0.0 R 0 = 0.0 , ỹ ∈ [a, b] and a, b ∈ R is the domain for the problem described in (1).The purpose of choosing the above set of ICs is to consider two types of population spatial distribution, i.e., homogeneous and heterogeneous environments to examine the role of diffusion in curtailing the infection in the aforementioned cases.The last two states, i.e., I A0 and R 0 are assumed to be zero because an epidemic usually starts with a relatively small number of initially affected people.Frequently, by the time the pandemic is identified, these people might not have reached the asymptomatic or recovered stage.
The initial conditions profiles are presented in Figs. 1 and 2. The illustration in Fig. 1 depicts the ICs (3), showcasing a uniform distribution of population across the domain for each of the sub-populations under investigation in this study.The susceptible people class is considered to be larger than the rest of sub-classes, while the asymptomatic and recovered population are considered to be zero.Figure 2 indicates the ICs (4) with exposed, susceptible, symptomatic, and individuals in super-spreading class concentrated around the center of the interval [−2, 2] and decreases exponentially towards the center (or origin) on both sides.The concentration of susceptible individuals at the origin significantly exceeds the concentrations of individuals in the exposed, symptomatic, and super-spreader classes.In addition, the concentrations of asymptomatic and recovered individuals are assumed to be zero.

Qualitative analysis of model
This section presents qualitative analysis of the reaction-diffusion COVID-19 compartmental epidemic model (1).We proceed to prove the basic mathematical properties of the model solution as follows.

Boundedness
One of the most important properties of an epidemic model is the solution boundedness.We take into consideration the approach described in 28 in order to analyze the solution boundedness of the problem (1).The result is given in the following Theorem.

Theorem 3.1
The solution of the model (1) i.e., (S(., t), E(., t), I 1 (., t), I 2 (., t), Proof In order to prove the desired result, add all equations in the model (1) Integrating over , using the well known Divergence theorem 29 and make use of no flux boundary conditions, which yield to

Invariant region
A positively invariant set for the system (1) is defined as follows:

Model equilibria and the basic reproductive number
For derivation of the threshold parameter known as the basic reproductive number, we used the well-known approached considered in 30 .Model (1) possess two equilibria i.e. the disease-free equilibrium (DFE) and the endemic equilibrium (EE) represented by 0 and Ŵ EE respectively such that: The EE is calculated as follows: with the following analytical values and Moreover, the basic reproductive number R 0 is computed as follows: The infectious classes in the proposed model (1) are E, I 1 , I 2 and I 3 , Henceforth, the vectors below present the transmission of newborn infections and the transitions between various classes.The Jacobian of above matrices are evaluated as: and N = � ζ in case of disease-free equilibrium.Thus, the associated next generation matrix is, Ŵ 0 = (S 0 , E 0 , I 1 0 , I 2 0 , I 3 0 , R 0 ) = (�/ζ , 0, 0, 0, 0, 0).Proof The Jacobian of (1) at the DFE point Ŵ 0 is and the characteristic polynomial associated to above matrix is given as follows: where, The values of c 0 , c 1 , c 2 in (6) are claimed to be positive under the condition R 0 < 1 .Further, c 1 c 2 − c 3 > 0 con- firming the Routh-Hurwitz conditions.Hence, the DFE stable locally when R 0 < 1 .

Local Stability of endemic equilibria
To discuss the local stability of endemic equilibria Ŵ EE , the model (1) is linearized at Ŵ EE = (S * , E * , I * 1 , I * 2 , I * 3 , R * ) .For this purpose assume that S( t, ỹ), Ē( t, ỹ), Ī1 ( t, ỹ), Ī2 ( t, ỹ), Ī3 ( t, ỹ) and R( t, ỹ) are minimal perturbation.The linearized form of the problem (1) is given by ( 8), and Given that the linearized system (8) has a solution in Fourier series form, then In above, k = nπ 2 , with n ∈ Z + , indicates wave-number for the node n.Using ( 9) in ( 8) yield to, www.nature.com/scientificreports/ The matrix V representing the variational matrix of ( 8) is where, The subsequent polynomial is, The relative coefficient values are; Vol The values of b i , where i = 1, . . ., 5 are given below: www.nature.com/scientificreports/The coefficients A i , i = 0, . . ., 4 , of P given by ( 13) are positive, if R 0 > 1 .Further, it also satisfies the Routh- Hurwitz stability conditions for a polynomial having degree five, i.e., Therefore, it is asserted that Ŵ EE is stable locally for R 0 > 1.

Global stability of the model
In the subsequent section, we will investigate the stability of the problem (1) in global case at the steady state Ŵ 0 using a nonlinear Lyapunov stability approach.
Theorem 3.3 If R 0 < 1 then the DFE Ŵ 0 of the system (1) is globally stable in .

The model's sensitivity analysis
Sensitivity analysis of epidemiological systems plays a crucial role in understanding the impact of system embedded parameters on disease incidence and prevalence.It helps to elucidate the significance of various parameters in shaping the dynamics of the epidemic.Due to the potential for errors in the collocation of data and uncertainties in parameter values, such analysis becomes essential for assessing the robustness of the respective model predictions.The sensitivity analysis provides insights into the reliability and stability of the model under different scenarios.To identify certain parameters that play a critical role and exert a significant impact on R 0 , computing their sensitivity indices proves to be effective.These parameters become prime targets for intervention strategies aimed at controlling the epidemic.The sensitivity index of R 0 relative t to the system parameter is defined as: where shows the system parameter involved in (1).Sensitivity indices of R 0 analytically represents as: Furthermore, to fulfill these analytical results, the sensitivity indices are evaluated numerically utilizing parameter values given in Table 1, and numerical indices are provided in Table 3.The above analysis indicates that the parameters β, r, β P , k 1 , k 2 , and ψ exhibit positive sensitivity indices, indicating their role in enhancing R 0 for larger values.Conversely, parameters ζ , ζ 1 , ζ 2 , η 1 , η 2 , and η 3 display negative indices, suggesting an inverse relationship with R 0 .Thus, an increase in these parameter values would decrease R 0 .Given their larger indices, these parameters, β and βP are the most sensitive.Small increments in these parameter's values could significantly elevate the value of R 0 .In biological context, the average number of contacts per person within a certain time interval is represented by the contact rate β .Elevated β values signify an increased probability of transmission during physical interaction.On the other hand, the high sensitivity of the parameter represents the transmission due super-spreader suggesting that the contribution of individuals with a high capacity for transmission, has a significant influence on the overall dynamics of the epidemic.

Numerical treatment: solution and simulation
The present part of the manuscript investigates the iterative solution of the system (1) using finite difference operator-splitting (FDOS) approximation technique presented in 17,18 .According to the Von Neumann stability criteria, the spatial and time step size for this purpose are considered as ỹ = 0.06 , and t = 0.02 days.The diffusivity coefficients have been taken as D 1 = 0.000050, D 2 = 0.00050, D 3 = 0.00010, D 4 = 0.0010, D 5 = 0.00010, D 6 = 0 .The proposed scheme steps are presented as fallow:

Solution scheme
The operator splitting iterative scheme is one of the efficient numerical techniques in the field of numerical analysis.This method is successfully applied for the solution of nonlinear partial differential equations in order to handle the complexity and non-linearity.This scheme is generally based on the splitting approach of differential operators into sub-operators.It results in the splitting of of problem under consideration into sub-problems corresponding to a particular physical phenomenon.In this investigation, the mentioned scheme is applied to solve the reaction-diffusion compartmental epidemic model for the dynamics of COVID-19 described in (1).Different population groups interact with one another and diffuse spatially in uni-direction requiring different temporal steps.Therefore, considering the operator-splitting approach in 17,18 , splitting the time dependent operator is useful which shift model (1) in to two sub-systems.The nonlinear reaction problem utilized for the time-step t 0 to 1 2 dt is Moreover, the diffusion problem in linear case utilized for time-step 1 2 dt to t n is as follows ( 14) www.nature.com/scientificreports/Now make use of finite-difference approximations, the time derivative with first order in ( 14) and ( 15) is approximated by as follows while the spatial derivative with second-order in the above system ( 15) is approximated by second-order central finite-difference described as follows: ξ stands for any of the variable S, E, I 1 , I 2 , I 3 , R. The iterative scheme for sub-systems ( 14) and ( 15) can be described as follows: and ( 15)

Simulation with discussions
This part presents simulation of the proposed spatio-temporal epidemic compartmental model (1) with the help of the numerical scheme discussed in ( 18) and ( 19) for the uniform and nonuniform initial conditions given by (3, 4).The impact of disease readmission rates β and β P has been depicted for different scenarios with and without diffusion.The time resulting simulation has been conducted for super-spreaders, exposed, symptomatic and asymptomatic individuals, under the initial conditions (3, 4) with as well as without diffusion at ỹ = 0.0 and ỹ = 1.0 .The choice of considering the spatial points ỹ = 0.0 and ỹ = 1.0 is according to the initial distribution of respective populations, which indicates the region having a high density of population.Moreover, the evolutionary trajectories are obtained for 700 days.The objective is to forecast future scenarios based on confirmed infected cases and to assess the effects of the aforementioned control interventions on exposed, super-spreader, symptomatic and asymptomatic individuals without as well as with diffusion cases.Only the ICs (4) have been considered for the system with diffusion, as ICs (3) assumes a uniform population distribution.A visual dynamics of the COVID transmission model ( 1) is presented to explore the effectiveness of diffusion in controlling the prevalence of COVID-19 infection utilizing the initial conditions (4).The dynamics of the infected human population without and with diffusion are presented in figure 3. It can be observed from these figures that curve peaks in all cases have a significant reduction in the presence of diffusion.Physically, it reveals that public gathering restriction plays a significant impact in minimizing the infection incidence.
Simulation for the initial conditions (3) at ỹ = 0.0 and ỹ = 1.0 We provide a visual depiction of the exposed, super-spreading, symptomatic, and asymptomatic population, both without and with diffusion, using uniform ICs (3). Figure 4 demonstrate the trajectories for the model (1) for ỹ = 0 and ỹ = 1 , utilizing the values outlined in the Table 2.A similar dynamics is noticed in each case as the IC (3) implies a uniform spatial distribution of the population.

Simulation based on ICs (4) at ỹ = 0.0 with different personal protection rates
This section accomplishes the impact of diffusion coupled with some of the model parameters for uniform and nonuniform ICs.The dynamics of the proposed model (1) for initial criterion given in ( 4) and under different interventions are illustrated in Figs. 5, 6, 7 and 8.These simulations are performed for both diffusive and non-diffusive scenarios at ỹ = 0.0 .Figures 5 and 6 depict the influence of β on the super-spreaders, exposed, symptomatic, and asymptomatic subgroups for both cases.Initially, the dynamics are examined for the tabulated value of β which is set at 0.5030.Variations in β correspond to changes in social contact intensity with increases and decreases representing relaxation and strengthening of social contacts, respectively.The impact of reducing β by 10% , 20% , and 30% is analyzed.It is observed that in the absence of diffusion, a 20% reduction in β results in a 53.01%decrease in infected individuals, while a 75.0%reduction is observed in the presence of diffusion.Furthermore, a 30.0%reduction in social contacts leads to a 75% decrease in infection without diffusion affect- ing exposed, super-spreaders, symptomatic, and asymptomatic infected individuals.Conversely, in the case of diffusion, a 93.01%reduction is calculated with a 30.0%decline in social contacts.Further analysis are presented in Tables 4 and 5. Consequently, from this analysis, it is evident that implementing isolation strategies in the presence of diffusion proves to be more advantageous and considerably contributes to mitigating the incidence of infection.Figures 7 and 8 present the impact of β P upon the population dynamics spatially at ỹ = 0.0 .Firstly, the visual dynamics are illustrated for the value of β P mentioned in Table 2 in diffusive as well as in non-diffusive models.The dynamics further analyzed for 10%, 20% and 30% reduction in β P .The simulation indicates that without diffusion case, the infected population in each compartment experiences reductions of 23% , 45% , and 64% , respectively.However, in the presence of diffusion, reductions of 35% , 64% , and 85% are observed in the corresponding compartments.The projected numbers are summarized in Tables 4 and 5.
� , � , � , This part of the paper presents the dynamical aspects of individuals in the aforementioned population groups at ỹ = 1.0 for 700 days.These results are graphically depicted in Figs. 9, 10, 11 and 12. Figures 9 and 10 describe the behavior of individuals in exposed, super-spreaders symptomatic and asymptomatic classes which are presented initially for β = 0.5030 without and with diffusion.To further assess the role on the respective infected classes, reductions of 10%, 20% , and 30% are implemented.In the absence of diffusion, as depicted by the initial concentration profile in Fig. 2, at ỹ = 1.0 , the low population concentration in compartments E, I 1 , I 2 , and I 3 results in a significant decline in the number of infected individuals with a 30% decrease in the effective trans- mission rate β .Conversely, when the population undergoes diffusion, the super-spreader, exposed, symptomatic and asymptomatic infectious individuals observed an increase for β = 0.5030 .However, the number of people in these classes effectively declined with a 30% decrease in β .Therefore, it can be concluded that implementing moderate social distancing policies is beneficial in reducing the number of infections in either scenario.The projected numbers for this case are presented in Tables 6 and 7.
Figures 11 and 12 visualize the influence of the infection transmission rate β P upon the dynamics super- spreaders, exposed, symptomatic, and asymptomatic infected individuals.The simulation are performed for varying β P , with a reduced by 10.0%, 20.0% , and 30.0%relative to the tabulated value.In case of no diffusion, the Vol:.(1234567890)   super-spreaders, exposed, symptomatic, and asymptomatic infected individuals experience a clear decrease due to the lower concentration level at ỹ = 1.0 , as depicted by the initial concentration profile in plot 2. Conversely, with the diffusion, the infectious individuals are increased at ỹ = 1.0 for β P = 0.7242 .With reductions in the aforementioned parameter, a decrease in infected individuals is observed in the respective compartments, with the lowest peak observed with a 30% decrease in β P .The corresponding plots agreed with the theoretical results, i.e., the solution stays positive and converges to the steady states throughout the domain.Further, the proposed numerical schemes preserve the positivity property.Moreover, the susceptible concentration is high at ỹ = 0.0 , according to initial profiles as given in Fig. 3. There- fore with baseline values of β and β P , the number of infected individuals in the respective compartments gets reduced with diffusion at ỹ = 0.0 and almost vanishes in the first 200 days.Thus diffusion will possibly curtail the infection in highly populated areas as it restricts public gatherings.

Conclusion
The proposed study is focused on the analysis of the dynamics of COVID-19 in a spatially heterogeneous case.The impact of some non-pharmaceutical interventions (such as personal protection, isolation, etc.) is observed with and without spatial effects.For this purpose, a spatio-temporal epidemic model is formulated, consisting of a system of partial differential equations describing the dynamics of populations with different disease statuses, Super-spreaders Individuals    in parallel with social distancing policy plays an important role in controlling and eradicating COVID-19 infection.
The present work can be extended to fractional diffusion problems using various operators for better understanding and control of the pandemic.

Figure 5 .
Figure 5. Simulation of exposed and symptomatic population under variation in β for diffusive and non diffusive case at ỹ = 0.

Figures 13 and 14
Figures 13 and 14 show the mesh plots of the reaction-diffusion model (1).The corresponding plots indicate the spatio-temporal evolution of the proposed model over the domain [a, b] × [0, T max ] , where a = −2 , b = 2 and T max = 700 days are taken for simulation purpose in order to investigate the long term behavior of the disease.The corresponding plots agreed with the theoretical results, i.e., the solution stays positive and converges to the steady states throughout the domain.Further, the proposed numerical schemes preserve the positivity property.Moreover, the susceptible concentration is high at ỹ = 0.0 , according to initial profiles as given in Fig.3.There- fore with baseline values of β and β P , the number of infected individuals in the respective compartments gets reduced with diffusion at ỹ = 0.0 and almost vanishes in the first 200 days.Thus diffusion will possibly curtail the infection in highly populated areas as it restricts public gatherings.

Figure 8 .
Figure 8. Dynamics of asymptomatic and super-spreader infected individuals under variation in β P for diffusive and non diffusive case at ỹ = 0.

Figure 9 .
Figure 9. Impact of parameter β over the solutions of exposed and symptomatic population with and without diffusion and ỹ = 1.

Figure 13 .
Figure 13.Mesh plots of susceptible, exposed, symptomatic and super-spreaders population in the presence of diffusion.

Figure 14 .
Figure 14.Mesh plots of asymptomatically infected and recovered individuals with diffusion.

Table 2 .
26ological details of the model parameters and respective numerical values.The values are taken from26.

Table 3 .
Sensitivity indices of model parameters versus R 0 .

Table 4 .
Projected outcomes of individuals in the respective infected classes of model (1) in non-defeasive case.

Table 5 .
Projected outcomes of individuals in the respective infected classes of model (1) in defeasive case.
Simulation of exposed and symptomatic infected individuals under variation in β for diffusive and non diffusive case at ỹ = 0.

Table 6 .
Projected peaks of the respective population classes of the model (1) at ỹ = 1.0 and without diffusion.

Table 7 .
Projected peaks of the respective population classes of the model (1) at ỹ = 1.0 with diffusion.